﻿using System;
using System.Collections.Generic;
using System.Data;
using System.Diagnostics;
using System.Linq;
using System.Numerics;
using MathNet.Spatial.Euclidean;
using MathNet.Spatial.Units;

namespace FDEE.Model
{

    /// <summary>
    /// 圆弧+直线+圆弧型预应力
    /// </summary>
    public class CSCTendonRecord : TendonBase
    {
        double AAnchorReserve; // A锚最小Reserve值
        /// <summary>
        /// 平直段长度
        /// </summary>
        double StraightLength;
        double CurveLength;
        double QZ;
        /// <summary>
        /// 曲线圆
        /// </summary>
        Circle2D Curve;
        Angle StAngle, EdAngle;
        /// <summary>
        /// 配置圆参数-三点定圆
        /// </summary>
        /// <param name="st"></param>
        /// <param name="ed"></param>
        /// <param name="md"></param>
        void SetCircle(Point2D st, Point2D ed, Point2D md)
        {

            var pt1 = st;
            var pt2 = ed;
            var pt3 = md;
            Point2D center;
            double radius;
            double x1 = pt1.X, x2 = pt2.X, x3 = pt3.X;
            double y1 = pt1.Y, y2 = pt2.Y, y3 = pt3.Y;
            double a = x1 - x2;
            double b = y1 - y2;
            double c = x1 - x3;
            double d = y1 - y3;
            double e = ((x1 * x1 - x2 * x2) + (y1 * y1 - y2 * y2)) / 2.0;
            double f = ((x1 * x1 - x3 * x3) + (y1 * y1 - y3 * y3)) / 2.0;
            double det = b * c - a * d;
            if (Math.Abs(det) < 1e-5)
            {
                throw new Exception();
            }
            double x0 = -(d * e - b * f) / det;
            double y0 = -(a * f - c * e) / det;
            radius = Math.Sqrt(Math.Pow(x1 - x0, 2) + Math.Pow(y1 - y0, 2));
            center = new Point2D(x0, y0);

            StAngle = Vector2D.XAxis.SignedAngleTo(st - center);
            EdAngle = Vector2D.XAxis.SignedAngleTo(ed - center);
            Curve = new Circle2D(center, radius);
            CurveLength = (EdAngle - StAngle).Radians * radius;
            QZ = ed.X;
        }
        /// <summary>
        /// 纵坐标
        /// </summary>
        double GetZ(double x_meter)
        {
            if (x_meter <= QZ)// 起点弧线
            {
                Point2D P = new Point2D(x_meter, Curve.Center.Y);
                double ll = Math.Sqrt(Math.Pow(Curve.Radius, 2) - Math.Pow(Curve.Center.X - x_meter, 2));
                return P.Y - ll;
            }
            else if (x_meter <= QZ + StraightLength * 0.001)
            {
                return GetZ(QZ);
            }
            else
            {
                return GetZ(QZ - (x_meter - (QZ + StraightLength * 0.001)));
            }
        }
        /// <summary>
        /// 摩擦损失
        /// </summary>
        /// <param name="x_meter"></param>
        /// <returns></returns>
        double GetFriction(double x_meter,out double ang_fromst)
        {
            double K = 0.0015;
            double mu = 0.25;
            double fpj = JackingForce * 1000 / Area;
            double xi = GetCurveL(x_meter);
            ang_fromst = GetAngleAbsChanged(x_meter).Radians;
            return fpj * (1 - Math.Exp(-(K * xi + mu * GetAngleAbsChanged(x_meter).Radians)));
        }
        /// <summary>
        /// 角度累计差
        /// </summary>
        /// <param name="x_meter"></param>
        /// <returns></returns>
        Angle GetAngleAbsChanged(double x_meter)
        {
            if (x_meter <= 0)
            {
                return Angle.FromDegrees(0);
            }
            else if (x_meter <= QZ)
            {
                return GetDirAngle(x_meter) - GetDirAngle(0);

            }
            else if (x_meter < QZ + StraightLength)
            {
                return GetAngleAbsChanged(QZ);
            }
            else if (x_meter < QZ * 2 + StraightLength)
            {
                var m1 = GetAngleAbsChanged(QZ);
                var m2 = GetAngleAbsChanged(QZ - (x_meter - QZ - StraightLength));
                return m1 - m2 + m1;
            }
            else
            {
                var m1 = GetAngleAbsChanged(QZ);
                return m1 * 2;
            }
        }
        /// <summary>
        /// 前进角度
        /// </summary>
        /// <param name="x_meter"></param>
        /// <returns></returns>
        public Angle GetDirAngle(double x_meter)
        {
            if (x_meter < 0)
            {
                return GetDirAngle(0);
            }
            else if (x_meter == 0)
            {
                return StAngle + Angle.FromDegrees(90.0);
            }
            else if (x_meter <= QZ)
            {
                var pt = new Point2D(x_meter, GetZ(x_meter));
                var vv = (pt - Curve.Center).Rotate(Angle.FromDegrees(90.0)).Normalize();
                return Vector2D.XAxis.SignedAngleTo(vv);
            }
            else if (x_meter < QZ + StraightLength)
            {
                return Angle.FromDegrees(0);
            }
            else if (x_meter < QZ * 2 + StraightLength)
            {
                return GetDirAngle(QZ - (x_meter - QZ - StraightLength)) * -1;
            }
            else
            {
                return GetDirAngle(QZ * 2 + StraightLength);
            }
        }
        /// <summary>
        /// 曲线长
        /// </summary>
        /// <param name="x_meter"></param>
        /// <returns></returns>
        public double GetCurveL(double x_meter)
        {
            if (x_meter <= 0)
            {
                return 0.0;
            }
            else if (x_meter <= QZ)
            {

                return Curve.Radius * GetAngleAbsChanged(x_meter).Radians;

            }
            else if (x_meter < QZ + StraightLength)
            {
                return GetCurveL(QZ) + x_meter - QZ;
            }
            else if (x_meter < QZ * 2 + StraightLength)
            {
                var m1 = GetCurveL(QZ);
                var m2 = GetCurveL(QZ - (x_meter - QZ - StraightLength));
                return m1 - m2 + m1 + StraightLength;
            }
            else
            {
                return GetCurveL(QZ) * 2 + StraightLength;
            }
        }

        //public Dictionary<string, double> ReportLine;


        /// <summary>
        /// 构造函数
        /// </summary>
        /// <param name="tMark">钢束编号</param>
        /// <param name="stage">张拉阶段</param>
        /// <param name="typeStrand">钢束类别</param>
        /// <param name="numStrand">钢绞线根数</param>
        /// <param name="isBothEnd"></param>
        /// <param name="JackForce">张拉力</param>
        /// <param name="straightLength_meter">直线段长度，单位m</param>
        /// <param name="St">三点定圆</param>
        /// <param name="Ed">三点定圆</param>
        /// <param name="Mid">三点定圆</param>
        public CSCTendonRecord(int tMark, int stage, eStrandClass typeStrand, int numStrand, eForceType isBothEnd, double JackForce,
            double straightLength_meter, Point2D St, Point2D Ed, Point2D Mid,double aAnchorReserve=0.75) :
            base(tMark, stage, typeStrand, numStrand, isBothEnd, JackForce)
        {
            StraightLength = straightLength_meter;
            AAnchorReserve = aAnchorReserve;
             ReportLine = new Dictionary<string, double>();
            SetCircle(St, Ed, Mid);
        }
        /// <summary>
        /// 实现长度算法
        /// </summary>
        public override void SetLength()
        {
            // 暂不考虑平曲线
            TheoricalLength = Math.Round(CurveLength * 2 + StraightLength, 3, MidpointRounding.AwayFromZero);
            //-------------------------------------------
         
            if (StrandsNum <= 18)
            {
                throw new NotImplementedException();
            }
            //else if (StrandsNum <= 19)
            //{
            //    AAnchorReserve = 0.8;
            //}
            else if (StrandsNum>19&&StrandsNum <= 22)
            {
                AAnchorReserve = 0.75;
            }
            else if (StrandsNum <= 27)
            {
                AAnchorReserve = 0.8;
            }
            else
            {
                throw new NotImplementedException();
            }

            double AAnchorB; // A锚最小B值
            if (StrandsNum <= 17)
            {
                throw new NotImplementedException();
            }           
            else if (StrandsNum <= 22)
            {
                AAnchorB = 0.325;
            }
            else if (StrandsNum <= 27)
            {
                AAnchorB = 0.43;
            }
            else
            {
                throw new NotImplementedException();
            }

            double AAnchorF; // A锚最小F值
            if (StrandsNum <= 17)
            {
                throw new NotImplementedException();
            }
            else if (StrandsNum <= 19)
            {
                AAnchorF = 0.075;
            }
            else if (StrandsNum <= 22)
            {
                AAnchorF = 0.08;
            }
            else if (StrandsNum <= 24)
            {
                AAnchorF = 0.082;
            }
            else if (StrandsNum <= 27)
            {
                AAnchorF = 0.085;
            }
            else
            {
                throw new NotImplementedException();
            }

            double PAnchorB; // P锚最小B值
            if (StrandsNum <= 17)
            {
                throw new NotImplementedException();
            }
            else if (StrandsNum <= 19)
            {
                PAnchorB = 0.72;
            }
            else if (StrandsNum <= 22)
            {
                PAnchorB = 0.9;
            }
            else if (StrandsNum <= 27)
            {
                PAnchorB = 1;
            }
            else
            {
                throw new NotImplementedException();
            }
            //-------------------------------------------
            if (StressType == eForceType.OneEnd)
            {
                // 管道长度=（理论长度）-（A锚B值）-（限位环厚度）-（P锚B值）
                DuctLength = Math.Round(TheoricalLength - AAnchorB - 0.02 - PAnchorB, 3, MidpointRounding.AwayFromZero);
                // 下料长度=0.2（工作长度）+（张拉长度）+（限位块厚度）+（锚头×N）+（理论长度）+（固定端预留）
                Length = Math.Round(0.2 + AAnchorReserve + 0.020 + AAnchorF * 3 + TheoricalLength + 0.2, 3, MidpointRounding.AwayFromZero);
            }
            else
            {
                // 管道长度=（理论长度）-2*（（A锚B值）-（限位环厚度））
                DuctLength = Math.Round(TheoricalLength - 2 * (AAnchorB + 0.02), 3, MidpointRounding.AwayFromZero);
                // 下料长度=2*(0.2（工作长度）+（张拉长度）+（限位块厚度）+（锚头×N）)+（理论长度）
                Length = Math.Round(2 * (0.2 + AAnchorReserve + 0.020 + AAnchorF * 3) + TheoricalLength, 3, MidpointRounding.AwayFromZero);
            }

            #region 导出计算过程--20211104

            ReportLine.Add("曲线长度", CurveLength);
            ReportLine.Add("平直段长度", StraightLength);
            ReportLine.Add("理论长度", TheoricalLength);
            ReportLine.Add("张拉端B值", AAnchorB);
            ReportLine.Add("限位环厚度", 0.02);
            ReportLine.Add("P锚B值", PAnchorB);
            ReportLine.Add("管道长度", DuctLength);
            ReportLine.Add("张拉端工作长度", 0.2);
            ReportLine.Add("张拉端预留长度", AAnchorReserve);
            ReportLine.Add("限位块厚度", 0.020);
            ReportLine.Add("锚头总厚度", AAnchorF * 3);
            ReportLine.Add("固定端预留长度", 0.2);
            ReportLine.Add("钢绞线下料长度", Length);
            ReportLine.Add("角度值", GetDirAngle(0).Degrees);
            #endregion
        }
        /// <summary>
        /// 实现伸长量算法
        /// </summary>
        /// <param name="dx"></param>
        public override void SetElongation(ref System.Data.DataTable ElongTable, double dx = 1e-3)
        {
            List<double> enlongation = new List<double>();
            double x0 = 0;
            double x1 = x0 + dx;
            double fpj = JackingForce * 1000 / Area;
            //double E = 196500;
            double E = 197000;
            double A = Area / StrandsNum;
            double dl;
            double rad_fromSt;
            while (x1 < QZ * 2 + StraightLength)
            {
                dl = GetCurveL(x1) - GetCurveL(x0);
                double xm = 0.5 * (x0 + x1);
                var f = fpj - GetFriction(xm, out rad_fromSt);
                enlongation.Add(f * A * dl * 1000 / E / A);

                #region ElongTable
                DataRow dr = ElongTable.NewRow();
                dr["No"] = string.Format("T{0}", MTendonMark);
                dr["Location"] = xm;
                dr["k"] = 0.0015;
                dr["μ"] = 0.25;
                dr["θst"] = rad_fromSt;
                dr["θed"] = 0.0;
                dr["fp0_st"] = fpj;
                dr["fp0_ed"] = 0;
                dr["fp_st"] = f;
                dr["fp_ed"] = 0;
                dr["fp"] = f;
                dr["dl"] = dl;
                dr["Ep"] = E;
                dr["△L"] = f * dl * 1000 / E ;
                ElongTable.Rows.Add(dr);
                #endregion

                x0 = x1;
                x1 += dx;
                if (x1 > QZ * 2 + StraightLength)
                {
                    x1 = QZ * 2 + StraightLength;
                }


            }
            Elongation = Math.Round(enlongation.Sum() * 1e-3, 3, MidpointRounding.AwayFromZero);

            #region 导出计算过程--20211104
            ReportLine.Add("张拉应力", fpj);
            ReportLine.Add("伸长量", Elongation);
            #endregion

        }

    }
}
